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Abstract 



We consider a nonlocal linear elastic wave model. We approximate 
the model using a spectral Galerkin method in space and analyze error 
! in such an approximation. We perform some numerical experiments 

^ ! to demonstrate the scheme. 

^ . 
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1 Introduction 

Integro-differential equations arise while modeling many problems in physics, 
biology and other areas which involve non-local diffusion/dispersal mecha- 
nisms j2j [lOl [H] . A nonlocal elastic wave model can be written as [IT] 



d 2 

—u(x,t) = P £u(x,t)+g(x,t) (t,x) enx(0,T), (1) 

where 

C(f)(x)= \ J°°(x -y)(j)(y)dy + a(x)(f)(x), 



n 



1 



(0,T), T > is the time interval under consideration, a(x) is a time- 
independent coefficient and g{x,t) is an inhomogeneity. Here we consider 
u(x, 0) = u , and ^m(x,0) — Vq, Q C K, where u(x,t) denotes the displace- 
ment field and d t u(x,t) represents velocity. Here the L 1 (M) kernel J°°(x) 
satisfies J°°(x) = J°°(—x) (symmetry). 

We assume that a(x) is a constant given by a = — j n J' OD (x)dx, and for 
simplicity we consider a normalized kernel function (J J°°(x)dx = 1). Then 

can be written as 

d 2 

—u(x,t) = pCu(x,t)+g(x,t) (2) 

where 

£<P(x)= [ J°°(x-y)(<f>(y)-<i>(x))dy. 
Jn 

Here we restrict ourself by J°° > (nonnegative) as well. 

Now let y — x = z, then dy = dz. Considering Q = M |3] in (T2]) we have 



I 1 = e J°°(x-y)(u(y,t)-u(x,t))dy = e / J°°{z) (u(x+ z,t) - u(x,t)) dz. 
Jr Jr 

Expanding u(x + z) in terms of Taylor series around z = we have 

h = £ J J°°(z) yu{x, t) + zu x (x, t) + ^—u xx (x, t) + - u(x, t)j dz. 

Thus we have 

Utt = eu x [ zj°°(z)dz + [ z 2 J°°(z)dz + + g(x, t). (3) 

As J(z) is symmetric, integrals with odd powers of z are zero. So takes the 
form u tt = C 2 u xx + C 4 u xxxx + ... + g(x, t). where C 2m = ^ f R z 2m J°°(z)dz; 
for each m — 1, 2, 3, - • • . Ignoring the higher order terms we have 

u tt = C 2 u xx + g(x,t). (4) 

The above equation is the familiar non-homogeneous wave equation, and is 
a first order approximation of the non-local equation (j2j). A comparison be- 
tween the operators u xx and Cu is well presented in [31 E] ■ Now following |3] 
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it is well understood that \\C\\ is bounded and C is negative semi-definite, if 
the kernel function J' 00 is considered non-negative, symmetric and normal- 
ized. These properties are important for the analysis we performed in this 
study 

Emmrich et. al. [IT] consider the non-local elastic model fl2]). They 
analyze mathematical and numerical solutions of the model. They compare 
the non-local model with that of a local continuum model. They claim two 
advantages of using this non-local model. They propose and implement a 
quadrature based numerical scheme considering infinite spatial domain. 

It is well understood from [U [§] that an integro-differential equation of 
type 02]) defined in the infinite domain can also be defined in a truncated finite 
domain [A, B] where A and B depend on the decay of the kernel function 
J°°(x). It is to note that a closed form formula to find suitable A and B is well 
presented in j9]. Thus the analysis and the approximation in a periodic Q = 
[0, 1] spatial domain can also be applied in any bounded interval [A, B] (and 
vice- versa). Considering the kernel of the convolution integral as gs{x) = 



y exp y— , the author in [9] formulate A = +y — 25 2 log(fev27r), 
and B = —A, where e > considered so that gs{x) > e. 

One can consider the model in a spatial periodic domain [3] as well. If we 
consider a spatially one-periodic initial function u(x,0), then for all x G R 
and t G M + u(x, t) = u(x + L,t). Then (j2j) can be written as 



u tt = (£u(-,t))(x)+g(x,t)= / J(x-y)(u(y,t)-u(x,t))dy + g(x,1$>) 



where J{x) = X^-oo J^^x—rL) for all x G [0, L], L > 0. We are interested 
to consider both periodic and non-periodic domain for spatial approximations 
of the model. 

Convolution models have been extensively investigated both theoretically 
and numerically [SI E], but a lot more yet to be done to get a highly ac- 
curate solver. Since an analytic solution can not be obtained always, an 
efficient and highly accurate numerical scheme has to be developed. Spatial 
approximation of this type of convolution models are interesting as well as 
challenging. The unknown function under integral sign and the nonlinear- 
ity involved in the model make the approximation more challenging. There 
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are various ways to handle such problems. In most cases scientists use some 
lower order schemes (with midpoint quadrature rules for integration [21 EI]) to 
serve their purpose. Thus there is much room for improvement and we find 
an interest of presenting and analyzing a higher order technique for space in- 
tegration. To the best of our knowledge, Legendre spectral Galerkin method 
has not been performed for this model problem yet and we find an interest 
to investigate such a scheme. 

In this article, we start with numerical accuracy analysis of a spectral 
galerkin approximation in section 2 followed by numerical implementations 
of the scheme in section 3. We present several computer generated solu- 
tions with relevent discussion in section 4. We implement the schemes in 
MATLAB. 

2 Accuracy of spectral Galerkin approxima- 
tion 

In this section, we look for spectral Galerkin scheme for space approximation 
and investigate convergence of such a scheme. Here we consider Legendre 
spectral polynomials. It is our goal to find u = u(x,t), such that 

uu(x, t) = Cu(x, t) + g(x, t),Wx G Q, 

and we seek for weak form to find u G L2(fl) such that 

(u u ,v) = (£u(x,t),v) + (g(x,t),v), V i)6l 2 (0). (6) 

We denote N be the set of all nonnegative integers. For any iV e N, P w 
denotes the set of all algebraic polynomials of degree at most N in Q. We 
define Legendre spectral Galerkin approximation of (T5]) as: Find 

u G Pjv, such that 

{v»,v) = {£u N (x,t),v) + {g(x,t),v), V v e¥ N (Q). (7) 
We define projection operator P^ such that 

(u,v) = (P N u,v), V v G Pat. 
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Then [T2J U\ 

\\u - P N u\\ < CN- m \u\ Hm(n) < CN- m \\u\\ Hm{n) 

where 



and 



and 



\u\n m (n) = I 



\ 1/2 

fci 



||W 

fc=minm,A r +l 

w|L i00 = max (Hit" 



{ll^lloo} , 

and C denotes a nonnegative constant which is independent of N. 

Let e N = u — u N be the error in spectral Galerkin approximation to the 
Legendre spectral Galerkin solution u N of (jBJ). Now from (E} and (J7|) we have 

((«-«%,«) = (£(«-«"),«). (8) 

Theorem 1. If u G is a solution of ffi), and u G S is a solution of (OJ) 
then the following inequality holds 

\\(u(;t)-U N (;t)) tt \\<CN- m \\u(.,t)\\m, 

for some C > 0. 

We need the following results to prove Theorem [TJ 

Theorem 2. [6] If J°°(x) > 0, J°°(-x) = J°°(x) for all x G Q and 

u(x) G M with u G L 2 (fl) then L is negative semi- definite. 

Theorem 3. f^lfj°°{x) > V x G R, J°°{x) G L 2 (R) and f R J°°(x)dx = 
1, then for any Q C M., L zs bounded and ||L|| < 2. 

Proof of Theorem LH We write 

e N = u — u N = u — Pnu + P/yu — = piv + 0jv. 



5 



Thus we rewrite flU]) as 

((9 N (; t)) u ,v) = (Cp N (;t),v) + (C9 N (;t),v). (9) 

Replacing v G by 9 N (-,t) G we get 

\((9 N {;t)) tt ,9 N {;t))\ < \{Cp N (;t),9 N (;t))\ + \(C9 N (;t),9 N (;t))\ 

and then applying theorem [2] (negative definiteness of C) 

\((9 N (;t)) tt ,9 N (;t))\ < \(C PN {-,t),9 N (;t))\. 

Thus applying theorem [3] 

W N (;t)) a \\\\9 N (;t)\\ <C\\ PN (;t)\\\\9 N (;t)\\ 

and then canceling common terms and bounds we get 

\\(0 N {;t)) tt \\ <CN- m \\u(;t)\\ Hm . 

□ 

3 Numerical implementation 

Here we implement spectral Galerkin schemes considering Legendre spec- 
tral polynomials. The Legendre polynomials satisfy three term recurrence 
relation: 

L (x) = 1, 
Li(x) = x, 

and 

(n + l)L n+1 (a;) = (2n + l)xL ri{x) - nL ri _i(x), n > 1, 

with 

1 2 
L k (x)Lj(x)dx = — — —Skj, 



2k + 1 

1 if k = j 
if k ^ j. 
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Here we consider Q = [— 1, 1] for numerical implementations and one may 
extend the scheme to any interval [-L, L], L > 0. We define 



N 



u N (x,t) = y~]a k (t)L k (x). 



k=0 

Replacing v G Pat by L k G Pat in §6§ we have 

(u%(x, t), L k {x)) = (Cu N {x, t), L k (x))+(g(x, t), L k (x)), V L k {x) G F N {Q). 

(10) 

Now 

N d 2 

(u%(x,t),L k (x)) = ^^aj(t)(Lj{x),L k {x)) 

k=0 

2 d 2 

~ 2kTidfi ak{t) > 

and 

N 

(g(x,t),L k (x)) = y]wig(xi)L k (xi). 

1=0 

where wi are Gauss weights and X\ are Gauss quadrature points. We define 

Cu(x, t) = Ciu(x, t) + C 2 u(x, t) 

= j ' J°°{x- y)u{y, t)dy - u{x, t) J J°°{x - y)dy. 

Now 

(jCu N (x,t),L k (x)) = (du N (x,t),L k (x)) + (£ 2 u N (x,t),L k (x)) 

where 

(Civ, N (x,t),L k {x)) = ^%(t) J J J°°(x - y)L k (x)L j (y)dydx 



and 



(C 2 u N (x,t),L k (x)) = J2 a j( t ) J L h {x)L 5 {x) J J°°(x-y)dydx. 
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Since we consider J°° such that 



I j°°(x) = 1, V iGfl, 
Jn 



we have 



(£ 2 u N (x,t),L k (x)) = Y]a.j(t) J L k (x)Lj(x)dx = 2fc ^ ^ k {t). 



Adding all these integrals f[TU|) can be written as 

2 d 



j2 n 



2k + l d^ ak ^ = ^2 a ^ I I J°°( x ~ y) L k(x)L j (y)dydx 



n Jn 



2k + 1 

for all k — 0, 1, • • • , iV. Thus (ITT]) can be written as 



OfeW + / g(x,t)L k (x)dx, (11) 



d 2 

M— a(t) = Aa(t) + 6(t), 

with a(0) = M _1 w , and o'(0) = M" 1 ^, A = Ai + A 2 and elements of At 
are defined from the continuous operator L\ as well as the elements of A 2 
are defined from the continuous operator of C 2 . Since M is a nonsingular 
diagonal matrix, the above second order ordinary differential equation can 
be written as 2 

^-a{t) = M- l Aa{t) + M~ l b(t). (12) 
dt 2 

We approximate the integrals by Gaussian quadrature formula 

/ / J°°(x - y)L k (x)Lj(y)dydx « S^wiLkfa) / JT ^ - y)Lj(y)dy 
Jn Jn l=0 Jn 

M M 

~ w i L k(xi)w m J OD (xi - X m )Lj{x m ), 



1=0 m=0 

and 

M 

g{x ) t)L k (x)dx m y^ j wig(xi,t)L j (xi), 
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where x\ are quadrature points, and w\ are quadrature weights. (fTTj) is a 
system of second order time dependent ordinary differential equations which 
can be solved using a simple finite difference scheme to approximate the 
coefficients a k (t). 

Depending on the nature of the kernel function, it is also wise to divide 
the integral domain into pieces. Then one can use small number of quadra- 
ture points to compute the integrals, since over each subinterval/subdomain 
J(x) will be flatter (if J(x) is a smooth function over Q.) We discuss the 
quadrature approximations in Section [A] 



4 Numerical experiments and discussion 

In this section, we exhibit several numerical results to illustrate the scheme. 
Here we use a simple scheme for time integration. We use MATLAB to serve 
our purpose. We approximate the time dependent system of equations ( 1T2|) 
using a central difference scheme 0. For simplicity we consider £1 = [—1,1] 
and the scheme can be extended to [-L, L], for any L > 0. We approximate 

m by 

n 3+l — On! -J- nJ- 1 

9 + = M~ l Aa> +l + M~ l b{tj) (13) 

at 2 

where a J = a(tj). This scheme is accurate pQ of order 0(At 2 + N~ m ). 

We know that p(A) < \\A\\ for any matrix A where p(A) stands for 
condition number of A. Also || f _ At a M -i A || < i-\\ d Ai-^A\y Now 

\A) k \ = | j J J°°(x - y)L k (x)Lj(y)dydx\ < max^jj J°°{x)\ J J L k (x)Lj(y)dydx\ 

<\jL k {x)dx\\jL J {y)dy\<AK, 



since 



t j( ,)dx- ' 2 3=0 



1 Mathematica and MATLAB builtin functions can also be used to solve the time de- 
pendent system of equations (|T2|) . 



9 



J°° is a normalized function and K = max^n \J' OD {x)\. We have 



\Ai\ 



N+l 



, E l^l 2 £ 



AT+1 



V 



16if = Vl6(iV+ l) 2 ^ 2 = 4(JV + 



and 



Thus 



I A, 



JV+l 



E I- 4 ?/ £ 2 



\ 



AT+1 

E 



2j - 1 



< 4. 



p|| < \\AxW + ||,4 2 || < 4 ((AT +l)K+l), 
and 1 < ||M|| < 4, shows that ||M _1 >1|| < 4 ((AT + l)K + 1) which is a 



bounded quantity for any finite N. So At 2 ||M 
Hence 

1 \ „ I 



l A\ 



P 



< 



< 



1 



I - At 2 M -1 v4 y - 11 1 - At 2 M~ 1 A " - 1 - || At 2 M-M| 
and thus by the definition of condition number of a matrix 

1 



as At 



< 1, 



0. 



P 



1, as At — >• and A" < oo. 



I -At 2 M~ 1 A / 

Thus the scheme ffl~3l) is unconditionally stable. 

In the following Figures [T]|2} we present the numerical solutions for var- 
ious choices of initial functions, p, and g(x,t) to demonstrate the scheme 
considering J°°(x) = <J M0 e - 400x ' 2 . Here we notice that the wave form is 
similar to the solutions of local linear wave model and the numerical solu- 
tions behaves well for any finite time preserving wave phenomena. Here 
in this study we have shown theoretically that the error in spectral approx- 
imation decays exponentially. More precisely, it is shown that if the kernel 
function and the solution underlying the model problem is smooth, then the 
errors obtained by the Legendre spectral scheme decays exponentially, which 
one desires from spectral method. 

There is a drawback of this approximation. Because of the presence of the 
convolution operator, the discrete operator is a dense matrix and as a result 
numerical computations become slower, needs huge storage for computations 
in multidimensional domain. 
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A Quadrature Approximation 



We define x k , k = 1, 2, • • • , M as quadrature points. Now in each point x h , 
we write <$fy as 

Qpu{x\t) = p J°°(x k -y)(u(y,t)-u(x k ,t))dy + g(x k ,t) 

which can be approximated by 
d 



32 M 



a(.r,t) = p^w^^ix 1 " -x i ){u{x\t)-u{x k ,t))+g{x k ,t),{U) 



in ! 

dt 2 



i=i 



where w l are quadrature weights. ( 1T4"|) is a system of second order ordinary 
differential equations, which can be solved by any standard techniques. 

One may also use a few quadrature points by subdividing the domain into 
several subdomains. We subdivide the domain fl into flj, j = 1, 2, • • • , N h 
subdomains such that Q = O n eac h °f the subdomains Qj, j = 

1, 2, ■ • • , Nh we define K points Xp k — 1, 2, • • • , K as quadrature points. 

Now in each point x k , we write ([2]) as 

d 



u(x k ,t) = p J°°(x k -y)(u(y,t)-u(x k ,t))dy + g(x k ,t) 
Jn 

m=l ^ rim 



which can be approximated by K point Gauss quadrature formula as 

B 2 N h K 

w u(x k ,t) = p^^<J«(x*-x{J(Ti( a 4,t)-t*(x*,t)) + ^ l ^) 

m=l i=l 

where are quadrature weights. (|T5|) can be solved using any linear ordi- 
nary differential equation solvers. 

Here, in Figure El we demonstrate a solution of the model in two space 
dimensions using mid-point quadrature formula for space integration followed 
by Euler's formula for time integration in a periodic (0, l) 2 space domain. 
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